Partitioning of energy in highly polydisperse granular gases 
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A highly polydisperse granular gas is modeled by a continuous distribution of particle sizes, a, giv- 
ing rise to a corresponding continuous temperature profile, T{a) , which we compute approximately, 
generalizing previous results for binary or multicomponent mixtures. If the system is driven, it 
evolves towards a stationary temperature profile, which is discussed for several driving mechanisms 
in dependence on the variance of the size distribution. For a uniform distribution of sizes, the sta- 
tionary temperature profile is nonuniform with either hot small particles (constant force driving) or 
hot large particles (constant velocity or constant energy driving). Polydispersity always gives rise 
to non-Gaussian velocity distributions. Depending on the driving mechanism the tails can be either 
overpopulated or underpopulated as compared to the molecular gas. The deviations are mainly due 
to small particles. In the case of free cooling the decay rate depends continuously on particle size, 
while all partial temperatures decay according to Half's law. The analytical results are supported 
by event driven simulations for a large, but discrete number of species. 

PACS numbers: 45.70.-n, 47.57.Gc, 47.45. Ab 



I. INTRODUCTION 

Granular media are an important and popular subject 
of current research which is owed partly to the striking 
phenomena they reveal and partly to their ubiquity in 
nature and in industry which makes a good understand- 
ing of their properties indispensable [U, [3, [1] . Of special 
interest are mixtures of different species, as real granular 
materials such as sand, gravel or seeds are rarely com- 
posed of identical particles. 

Starting with Jenkins and Mancini 0, Q binary mix- 
tures and in particular their kinetic temperature and 
transport properties received considerable interest @, 0, 
i, i, M, [U, m, El, [II- These studies confirmed that 
equipartition of energy is indeed violated in granular bi- 
nary mixtures, an observation that was first made in ex- 
periments by Losert et. al. p^ . Polydisperse granular 
mixtures, i.e., mixtures composed of more than tw o ty pes 
of particles were studied much less [H, [13, [H, 0, [lo, lHI 
although they are closer to realistic systems. In partic- 
ular, Dahl et. al. [l3l and Zhi-Yuan et. al. [2lj sim- 
ulated mixtures of particles with a distribution of sizes, 
and Lambiotte et. al. [l^ discuss mixtures of Maxwell 
molecules with varying coefficients of restitution. 

Out of the many fascinating phenomena inherent to 
granular mixtures and the observables that are necessary 
to understand them, we will focus on the partitioning of 
energy and how it evolves in time, both in the homoge- 
neous cooling state (HCS) and in homogeneously driven 
systems. Even though, in this paper we will first develop 
the machinery to deal with an arbitrary number, X, of 
species, we will eventually go one step further and con- 
sider highly polydisperse systems, where no two particles 
are alike but instead possess properties that are drawn 



from continuous probability distributions. 

In the following three sections we give a short intro- 
duction to the model and methods we use. In section fVl 
we investigate the temperature in a highly polydisperse 
system, characterized by a continuous distribution of par- 
ticle sizes. We finish with a brief conclusion and delegate 
all technical material to the appendices. 

II. MODEL AND OBSERVABLES 

In order to model a polydisperse granular gas, we con- 
sider mixtures of X different species of smooth inelastic 
hard spheres. Each species a = 1, 2, . . . , A" consists of 
Na oo identical particles, such that the concentra- 
tions Xa ■= Na/N {N = J2a ^a) as wcll as the density 
n = N/V remain finite as — ^ oo. Collisions between 
particles are assumed to be instantaneous and the parti- 
cles move freely between collisions. Because of the van- 
ishing collision time collisions of more than two particles 
can be neglected, i.e. the dynamics is determined by 
two particle collisions. The inelasticity is described by 
a velocity independent coefficient of normal restitution, 
eQ/3 G [Oil]; which may depend on the pair of species 
a, /3 = 1, 2, . . . , AT that the colliding particles belong to: 

h ■ ^ -eapfl- Vi2, (1) 

where Vi2 = Vi — V2 is the relative velocity of the collid- 
ing particles at contact before the collision and the 
corresponding quantity after the collision. The unit vec- 
tor h points from the center of particle 1 to the center of 
particle 2. Apart from the mutual coefficient of restitu- 
tion Eap, the species may also differ in mass TOq and in 
size (radius) Oc,. 
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The collision law [cq. ((T])] together with conservation 
of momentum determines the postcollisional velocities v[ 
and v'2 uniquely in terms of the precollisional ones {vi, 

V2)- 



Vi= Vi- 



m2 



7711 + "I2 

mi 



mi + m2 



(1 + ei2)(n • vi2)n, 
(1 + ei2)(n ■ vi2)n 



(2) 



As we consider smooth spheres the tangential component 
of the relative velocity {vi2 x n) remains unaffected. 

Due to the inelasticity, the particles suffer an energy 
loss during collision, i.e. the gas will cool down. To com- 
pensate for this energy loss, one can provide the system 
with external energy. We will restrict ourselves to volume 
driving [2^ : With a given frequency /dr random kicks 



Pt Pi +PdAtit) 



(3) 



are applied to all particles individually {pt = rriiVi). The 
strength of the kicks is controlled by pdr while the com- 
ponents of are drawn form a white noise source: = 

and £.f{t)^'^{t') = S^jS^^Sit - t'). The time between two 
driving events is taken to be small compared to the time 
scale on which the gas would cool down without energy 
supply. 

When considering X-component mixtures, the driving 
strength pdr may in general be a function of the particle 
species pdr = Pdf There are several experimental meth- 
ods (both m. D ~ 2 and D = 2>) that one can hope to 
describe approximately by volume driving: Shaking on a 
rough plate [2^ , electrostatic |23. 25| or magnetic [29, 
excitation, fluidisation by air 27, 23| or water |29l |. As 
it is not obvious how to best describe the driving of all 
these experiments theoretically, we propose the following 
three simple mechanisms: 

i. force controlled driving, assuming that all particles 
experience the same force (p^r = Pdi)i 

ii. velocity controlled driving, assuming that all par- 
ticles get velocity kicks of the same magnitude 
(p^^ oc m„) and 

iii. energy controlled driving, supplying every species 
on average with the same energy (p^^. (x m]/^). 

The first two mechanisms combined with an additional 
viscous drag force oc rjv are also discussed in the context 
of binary mixtures by Pagnani et. al. [13 . Our hope IS 
that the results discussed below may help to clarify the 
experimental conditions. 

The basic quantity of interest is the one-particle ve- 
locity distribution, fa{v)d^v, of species a which is re- 
lated to the one-particle distribution fa{r,v)d^rd^v by 
fa{v) = J fa{r,v)d^r. As an example, consider species 
that differ in mass, so that the one-particle velocity dis- 
tribution is explicitly given by 



N 



.{5{v-v,))d^v, 



where the angular brackets (•) denote the average over 
the A^-particle distribution function. It is normalized 
such that 

J d^'vUiv) ^ 7V„ and J ^""^U-") = 

a 

The partial granular temperature for species a in D space 
dimensions is defined by 



Jd^vUv) ■ 
(4) 

The mean temperature, T — J2a ^aTa, is then just given 
by the mean kinetic energy 

= 1 y / d^« uv)"^ = 1 y ^K^). 

a •' i 

The above definitions are easily generalized to other 
species characteristics, e.g. different size or different co- 
efficients of restitution: The indicator function, S„i.^„i^, 
just has to be replaced by the corresponding one. 

Our main emphasis in this paper are particles whose 
properties depend on a continuous variable a e R that 
follows a prescribed probability distribution dfi(a), i.e. 



E 



N 



dax{a) = / d/i(a). 



The temperature becomes a continuous function Ta 
T{a) whose mean and variance is given by 



T = / T{a)d^i{a) , AT = - T 
with T^ ^ I T^{a)dfi{a 



(5) 



In our example of a distribution of masses, a = m the 
one-particle velocity distribution, f{m,v)d^vd'm, is de- 
fined by 

N 

f{m,v) = ^(5(mj - m){5{v - Vi)). 



III. ANALYTICAL THEORY 

The time evolution of the temperatures is computed 
with the help of the pseudo Liouville operator formalism. 
For details sec e.g. rcfs. [s^ and [3l|. In this framework 
the time evolution of an observable A is given by the 
equation 



dt 



(A) = {^LA), 



where iH denotes the pseudo Liouville operator. 
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The pseudo Liouville operator for the driven hard 
sphere gas consists of three terms. The term jLq de- 
scribes free streaming which docs not affect the temper- 
ature, the term iLh accounts for driving and i-C/ for 
interactions between particles. In a gas consisting of X 
different species one obtains 



H 



X a 

EE 

a=l fj=l 



where i-Cq/j accounts for interactions between particles of 
species a with particles of species (3. For the evolution of 
the temperature of a particular species, only interactions 
with participation of that species play a role; collisions 
between particles of other species do not have a direct in- 
fluence. Given a discrete number, X, of different species, 
the temperature of species a, eq. (g]), develops in the 
following way 



D d 
2"di 



Ta = {iLHEkm{a)) + ^(i'CQ/3i?fc,„(a)). (6a) 

/3=1 



Given a continuous distribution d^{a) of a parameter a, 
one obtains 



D d , , 



{iLnEkinia:)) + J {iLaf3Ekin{a))d^{l3). 

(6b) 

At this point we would like to stress that the above equa- 
tions hold subject to arbitrary initial conditions Ta{t = 
0). The a priori assumption of a (quasi-)stationary state 
that is required for some of the hydrodynamic theories is 
not needed here. 

For a hard core potential the interaction terms iHap 
separate into a sum of two particle interaction opera- 
tors i£jaf3 = \ 'Ylik I ^"^^ap "^ith ouc particlc belonging to 



by 



dr 



,)\{d^rkA''vk---\{A^r^A^r 



k=l 



with rij the distance between particles i an j. 

We assume that the particles are uniformly distributed 
in space, that the species are well mixed and that velocity 
correlations between different particles can be neglected 
(molecular chaos assumption) . Under these premises the 
A^-particle distribution function fN{{ri},{vi},t) factor- 
izes into a product of N single particle distribution func- 
tions f(r,v,t). In a monodisperse system, the single par- 
ticle distribution hmction can be written in rescaled form 



T{t)D/2 



both in the homogeneous cooling state as well as in the 
stationary state of a driven system (3^ . 

Contrary to elastic gases, a Gaussian distribution is 
only an approximate solution in the inelastic case. Devi- 
ations have been studied extensively for driven and un- 
drivcn monodisperse systems. Investigations have shown 
that while a Gaussian approximation is quite good in the 
range of typical velocities, hi gh v elocities are overrepre- 
sented in granular gases [3^ l33l. [33 . Issl \3^ . In ref. 
qualitatively similar deviations have been found for freely 
cooling binary mixtures. The corrections have, however, 
only little influence on the temperature and the cooling 
rate |32l |. Thus, we make a Gaussian ansatz for the ve- 
locity distribution of a single species a with temperature 
Tq,. The iV-particle distribution for a mixture with X 
components then follows: 



species a, the other one to species /?. For the operator fN{{fi} {vi} t) cx ]^e~ 



i7^''9 one obtains 

otfj 

Jkl) 



2Ti(t) 



i=l 



Nx 

n 



(7) 



ilap ■■= -{vki ■ n)e{-Vki ■ n)S{rki - ak - ai){b^^j^ - 1), 

where b'^^j^ is the operator replacing the particles' veloc- 
ities before collision by their values afterwards according 
to equation 

When calculating the phase space average, one has to 
take into account the excluded volume effect which arises 
due to the fact that particles cannot overlap. Conse- 
quently, the phase space element in D dimensions is given 



In undriven systems, the HCS is maintained only for a 
certain time until velocity correlations develop and clus- 
ters form because of the system's instability against den- 
sity fluctuations (37l . [38l . [soj . In inelastic mixtures clus- 
ter formation is additionally accompanied by the onset 
of segregation S^l . Therefore our results will in this 
case be limited to the initial development. 

Using the distribution function, cq. ([7]), evaluation of 
the term {iLapEkinicn)) yields (cf. appendix [A]) : 



{i^apEkin{oi)) = —2xi3fJ.apGapi 



I 



(,J/3 — -La, 

rria mp 



(8) 



with the reduced mass '■= i^ai^p/i'm-a + i^p)- The other constants are given by 



A R ._l(l + e^ 
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Gap ■■= 4(aa + af3)n. Xq/3 for £» = 2, 



and 



Gap ■■= 8{aa + ap) n. Xap for D = 3 

V ITla 

where Xap is the value of the pair correlation function 
gapif) at contact. In the following, we will use the ap- 
proximation Xra = 1 which is well justified for dilute 
systems. 

The terms in equation ([8]) have a direct physical in- 
terpretation: The factor before the square brackets de- 
fines an effective collision frequency ujap of particles com- 
ing from possibly different species with different tem- 
peratures. The first term inside the brackets accounts 
for the dissipation in collisions between a and f3 parti- 
cles while the second term describes the heat flux be- 
tween species with different temperatures which tends to 
equalize the two temperatures. This term is the only 
one present in mixtures of elastically colliding particles, 
where it ensures equipartition. The difference to the elas- 
tic cases consists in the dissipative terms. As the cooling 
rates oc GapAap are in general different for each species 
and are completely independent from the rate of energy 
exchange oc Gap Bap they constantly drive the system 
away from equipartition. The new quasi-stationary state 
is then no longer characterized by equipartition but by 
equal cooling rates Ta/Ta = Tp/Tp A related in- 
terpretation has been given before by Alam and Luding 
|13| . Moreover it is also apparent that driving the system 
will in general not be sufficient to restore equipartition 
as was first shown by Barrat and Trizac [8|. For the 
special case of an undriven system that already reached 
its quasi-stationary state, equation (j6ap is equivalent to 
equation (2.4) of ref. l20la. 

The driving power Ha which was formally written as 
Ha — {i^HEkinioi)) in equations (|6a[l and \&o\ can 
be more easily calculated directly form the definition 



[eq. @] 



Ha 



fM,y/2ma. 



In particular we get for (i) force controlled driving H^^ oc 
l/rUa, (ii) velocity controlled driving H^ oc ma, and (iii) 
constant energy input H'^ independent of nia ■ 



IV. SIMULATIONS 

In order to test our analytical theory we performed 
complementary computer simulations based on an event- 
driven (ED) algorithm [Jlj. Although our code can easily 
handle up to 10^ particles, we usually found 10^ particles 
per species sufficient for the measurements reported here. 
Because of the extremely low densities used in this paper, 
we hardly ever need to take care of the inelastic collapse 
occurring in ED-simulations. If necessary we use the 
method of ref. to avoid inelastic collapse. 



For monodisperse systems, the minimal cluster size Lc 
can be derived from a hydrodynamic stability analysis 
m, HI]. To keep our systems from clustering, we chose 
a system size L < Lc/6. Although Lc will certainly be 
somewhat different for polydisperse systems, we found 
no indications for clustering or segregation in our simu- 
lations. 

As mentioned above, our simulations include volume 
driven systems. In this context it is necessary that the 
simulation process takes the conservation of momentum 
into account. To do so, a driving event always concerns 
two particles at the same time |53| . One of these parti- 
cles, say particle 1, is chosen at random. The neighbor- 
hood of this particle is examined to find the particle, i, 
closest to the first one. Particles 1 and i are then kicked 
at the same time t. While a momentum increment Pdr^(t) 
[see eq. ([3])] is added to particle 1, it is subtracted from 
particle i, i.e. 



Pi 
Pi 



Pi 

Pi 



of 



In that way momentum is conserved on length scales 
a mean particle separation, i.e., £ oc n^^/^ . 

The simulations were performed in two steps. Initially 
the particles were placed on a grid and random veloci- 
ties drawn from a Gaussian distribution were assigned to 
the particles. In the first half of the simulation, all co- 
efficients of restitution were set to unity and the elastic 
mixture was simulated for about 120 collisions per par- 
ticle to generate a well mixed state. In the next step the 
desired inelasticities were switched on and the tempera- 
tures were recorded until reliable estimates for the sta- 
tionary values of the observables could be obtained. For 
the driven systems we chose the driving frequency /dr to 
be approximately the same as the collision frequency at 
the desired stationary temperature Too. As a compro- 
mise between computational efficiency and the desire to 
reduce temperature fluctuations due to rare but strong 
driving events this choice of driving frequency was also 
found satisfactory by Bizon et. al. |44l |. 



V. HIGHLY POLYDISPERSE SYSTEMS 

Many real granular systems are highly polydisperse 
with no single particle being identical in shape and size 
to another one. To account for a high degree of poly- 
dispersity we generalize the considerations for polydis- 
perse mixtures to mixtures containing "inflnitely" many 
species. In principle, a variety of scenarios can be thought 
of and treated within our analytical approach. Here, we 
will restrict ourselves to the relatively simple case where 
the particles' radius is uniformly distributed in a range 
[i?i,i?2]; the particles all have the same mass density p 
and all restitution coefficients are equal eap = £■ We 
furthermore choose units such that p = I. 

The following questions are of particular interest. Is 
there a stationary temperature profile, T(a), if the sys- 
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tern is driven? If so, how docs this function reflect the where the nonhncar integral operator F is given (in D = 
properties of the distribution of radii? How does the fore- 3) by 
ing mechanism affect the stationary temperature profile? 
How does the system cool freely if undriven? 

Combining equations (j6bp and ([8]) leads to the follow- 
ing integro-differential equation for the temperature of 
species with radius a: 

^±Tia)=Hia) + F[T]{a) (9) 



F[T](a) 



l^/6 



i?2 — ^1 



drx 



r3(a + r)2 jf(a) T{r) 



ra ^2 



[T{r)-T{c 



When the system is driven constantly in time, we ex- 
pect a stationary temperature profile Too (a) = T{a,t — > 
oo), to develop. If this is correct, it should be given as 
the asymptotic solution of equation © with the left hand 
side set to zero: 



f[T^](a)^-H{a) 



(10) 



In general T^c depends not only on a but also on the two 
parameters i?i,i?2 of the distribution of radii. By scal- 
ing all radii with Ri, one observes that (up to a scale 
factor) Too depends only on the ratios a* = a/Ri and 
R = R2/R1, but not on the absolute values. Alterna- 
tively, we choose a* and the relative width of the dis- 
tribution A = 2(i?2 ~ Ri)/{R2 + Ri) as independent 
variables: T^o = Too (a*, A). 

We solved the above nonlinear integral equation (fTUl) 
numerically by applying Banach's fixed point iteration 
(for details see appendix[Bj. We always found a solution, 
confirming that a stationary temperature profile is indeed 
reached for asymptotically long times. 

Independently, we performed event driven simulations 
and measured all the partial temperatures T{a,t). The 
amount of simulation time needed for sufficiently good 
statistics quickly rises with the number of species. To 
this end, we checked if a highly poly disperse system can 
be approximated by a polydisperse mixture with many 
species such that there is still a considerable number of 
particles for each species. Considering equation ((6a|) for 
increasing numbers of species we found that the temper- 
atures considered in this paper rapidly converge. Fig- 
ure [Ua) shows how mixtures of respectively three and 
five species compare to the result for a continuous distri- 
bution. From these results we conclude that considering 
X = 20-30 species for the simulations should yield results 
practically indistinguishable from the highly polydisperse 
case. 

In FIG. m we show the stationary temperature 
Too (a*, A) as a function of particle radius a* for the 
three driving mechanisms proposed in section |TT1 The 
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FIG. 1: a) The stationary temperatures [eqs. (|6a|) & ()6b|l ] 
in a three component (open disks) and five component (filled 
disks) mixtures compared to those of a highly polydisperse 
mixture for energy controlled driving H"'^ = 10""^ at density 
71 = 5 X 10~* and coefficient of restitution e = 0.9. b) Inverse 
cooling time uiq in a two dimensional system for a uniform size 
distribution of width R = 3, coefficient of restitution e — 0.9 
and density n = 2 x 10"''. The symbols denote simulation 
results for X = 30 species each with 10'* particles, while the 
solid line is the solution of eq. (|12p . 



rough trends can be understood from the following qual- 
itative arguments. Force controlled driving H^'^{a*) oc 
l/m(a*) oc a*~^ is dominant for small particles so that 
one expects the partial temperatures, Too(a*, A), to de- 
crease with increasing size a* . This is indeed born out by 
the solution of the integral equation (fTU]) and supported 
by simulations, which are seen to agree well with the the- 
oretical result. Velocity controlled driving H'^'^ oc m{a*) 
is dominant for large particles so that we expect the par- 
tial temperatures to increase with increasing size of the 
particles, as is indeed observed in FIG.[2] Finally, for the 
energy controlled mechanism, H°'^{a*) = H, is indepen- 
dent of the particle size, nevertheless Too (a*, A) depends 
weakly on a* . One has to keep in mind that all the species 
interact and that this will lead to nontrivial conditions 
of stationarity as in the binary case. These effects are 
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FIG. 2: (Color online) The stationary temperature in a three 
dimensional driven system for a size distribution with R2 = 
3i?i, a coefficient of restitution e = 0.9 at a density n = 
2x10"'*. Force controlled driving i?'^'^ (a) = 1.875x 10"^/m(a) 
(solid, red), energy controlled driving H'"^ = 1.875 x 10^^ 
(long dashed, blue) and velocity controlled driving H^'^{a) — 
1.875 X 10~^m{a) (short dashed, green). For the simulation 
data (symbols) a system of 20 different particle species with 
10'' particles each was used. 



responsible for the precise functional form of the temper- 
ature profile which goes beyond the simple rough trend 
for all three driving mechanisms. The same trends for 
force controlled versus velocity controlled driving have 
been found by Pagnani et. al. [l^ in the case of binary 
mixtures. 

Abate and Durian [2^ discuss several systems that, al- 
though they are comprised of only two to five particles 
come close to our definition of highly polydisperse sys- 
tems in that no two particles are alike. Two spheres of 
different sizes show a marked increase in the tempera- 
ture ratio with increasing size ratio. This would roughly 
correspond to our results for velocity controlled driving 
but the authors of ref. [H observed a complicated two 
particle interaction. Moreover, they considered a system 
with five different spheres of the same size but different 
densities. Based on the results from binary mixtures one 
infers that the effects of different masses is much stronger 
than that of different sizes. If this reasoning is valid the 
weak dependence of the temperature on the mass would 
correpond to energy controlled driving in the present pa- 
per. 

Within our approximation scheme, the partial temper- 
atures (i.e., the temperature profile). Too (a*), determine 
the one-particle velocity distribution according to 



f{a,v) = 



N 



i?2 ~ ^1 



m(a) 



D/2 



The total velocity distribution, f{v)(Pv is thus given by 



R2 



f{v) = / daf{a,v). 



(11) 



This function is in general not Gaussian, not even for 
an elastic molecular gas with many different species. In 
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FIG. 3; (Color online) Stationary velocity distribution 
[eq. (|lip ] in a three dimensional driven system; parameters 
and symbols as in FIG. (2] for comparison the velocity dis- 
tribution of an elastically colliding, molecular gas (dashed- 
dotted, black) and a gaussian fit to the central part of the 
distribution (thin dotted line) are also shown. 
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FIG. 4: Stationary velocity distribution relative to the elastic 
gas and separately for the two halves of smaller (full line) 
and larger (dashed line) particles; parameters as in FIG. [2] 
a) force controlled driving, b) energy controlled driving. 



FIG. [3] we show the total velocity distribution as given 
by equation PT]) . The elastic system (dashed-dotted) is 
compared to the inelastic gas with different driving mech- 
anisms. In comparison to the molecular gas the tails of 
the velocity distribution can either be overpopulated, as 
observed for force controlled driving (solid line), or un- 
derpopulated for energy (long dashed) or velocity con- 
trolled (short dashed) driving. 

To clearly see the difference to the elastic case, we plot 
in FIG. Slthe velocity distribution relative to the elas- 
tic gas. We furthermore separate the particles into two 
halves, one with the smaller and one with the larger par- 
ticles. The strongest deviations are clearly in the tails 
and solely due to the small particles. The velocity distri- 
bution of the large particles has almost the same form as 
in the elastic gas, except for very small velocities. Force 
and energy controlled driving are almost mirror images 
of each other — even for the detailed structures at small 
velocities. 

How does the temperature profile. Too (a*. A), reflect 
the prescribed distribution of radii? The latter is char- 
acterized by a single parameter, the relative width A, 
which can take values < A < 2. In FIG. Owe show 
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FIG. 5: (Color online) a) The mean temperature T [eq. ([SJ] 
and b) the relative temperature variance AT^ = T'^ /T — 1 
as a function of the relative width of the size distribution 
A for the three driving mechanisms H^'^{a) = 10~^/m{a) 
(solid, red), H'^'ia) = 10"^ (long dashed, blue) and H'"'{a) = 
10~^m{a) (short dashed, green), coefflcient of restitution 
e — 0.9 and density n = 5 x 10~*. The mean temperatures in 
(a) are rescaled such that they agree for a relative width of 
A = 1. 

the mean temperature T and the temperature variance 
AT2 ■.= T^-T^ (see eq. ©) as a function_of In 
FIG. [5fa) we show the mean temperature T(A)/r(l), 
scaled such that they coincide at A = 1 . Surprisingly the 
dependence is nonmomotonic for different driving mech- 
anisms: whereas T(A) increases with A for force and 
velocity controlled driving, T(A) decreases with A for 
energy controlled driving. The strongest variation is ob- 
served for force controlled driving. The corresponding 
variance of the temperature profile [FIG.Efb)] increases 
trivially with A. The variance for velocity controlled 
driving is almost an order of magnitude larger than for 
the other two driving mechanisms. 

We next consider the freely cooling case {H{a) = 0) 
[mil. We expect Haff's law fis^ to hold also for T(r, t) 
and hence make the ansatz 

T{a, t) oc LJo{a)^'^t~'^ 

for large times. This leads to an integral equation for the 
inverse cooling time 

- Du;^\a) ^ F[u;^']{a) (12) 

Similarly to T^o, the decay rate too depends only on 
a* ~ a/Ri and R = R2/R1, or alternatively A, but 
not on the absolute values: loq = LL!o(a*^A). The 
above integral equation is solved numerically by subse- 
quently applying Banach's fixed point iteration and New- 
ton's method (for details see appendix [B|) . To extract 
ujQ^a* , A) from the simulations, we performed simula- 
tions with X = 30 species, measured all partial tem- 
peratures and fitted them to Half's law. The resulting 



decay rates are plotted in FIG. [T]Jb). The rate is seen 
to be a monotonically decreasing function of a* , however 
the dependence is weak. Since the coefficient of restitu- 
tion is the same for all particles, this is a pure size effect, 
implying that smaller particles relax faster than larger 
ones. The simulation data are seen to agree well with 
the theoretical results, but show a considerable scatter. 
This is most likely due to the difficulty in fitting the data 
to Haff's law, given the uncertainty in time scale, when 
the asymptotic decay applies. 



VI. CONCLUSION 

We examined the partitioning of energy in highly poly- 
disperse mixtures of smooth hard spheres. The proper- 
ties of the particles, such as mass, radius or coefficient 
of restitution, are chosen from a continuous distribution 
giving rise to a corresponding continuous temperature 
profile. The latter has been computed approximately, 
generalizing previous approaches of mixtures with sev- 
eral species. The analytical theory leads to a nonlinear 
integro-differential equation for the time dependent tem- 
perature profile, which has been solved numerically. 

Our results are supported by event driven simulations 
for mixtures with X = 20 — 30 species. The good 
agreement between ED simulations and the analytical 
theory indicates that the assumptions of homogeneity 
and molecular chaos that are fundamental to the the- 
ory are indeed observed in the simulated system. The 
direct simulation monte carlo (DSMC) method [4^, oth- 
erwise well suited for dilute (granular) gases (see e.g. 
113, S, Si, HI), would not have been able to show this 
as it ensures both homogeneity and molecular chaos by 
construction. 

As a specific example we have studied a uniform size 
distribution in detail. We showed that a highly polydis- 
perse mixture still obeys Haff's law during free cooling. 
The distribution of sizes gives rise to a nonuniform dis- 
tribution of cooling rates, such that the smaller particles 
are cooling faster. 

A driven system relaxes to a stationary temperature 
profile which is in general nonuniform. Depending on 
the driving mechanism, its weight can be predominantly 
at small or large particles. If the particles are driven 
by a constant force, then the smaller particles are hot- 
ter. If the driving process supplies either a constant en- 
ergy or velocity, then the larger particles have a higher 
temperature. The temperature profile reflects the distri- 
bution of radii, characterized by the relative width A. 
The variance of the temperature increases with A, as 
one would expect, whereas the mean temperature can ei- 
ther increase (constant force driving) or decrease with A 
(constant energy supply). 

This strong dependence on the driving mechanism is 
also observed in the velocity distributions. For a poly- 
dispcrse system, these are in general weighted sums of 
all partial distributions and hence in general not Gaus- 
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sian, even if the partial distributions are Gaussian like in 
an elastic gas. The velocity distribution in an inelastic, 
driven gas can have either overpopulated or underpopu- 
lated tails, as compared to the molecular gas. Further- 
more, the effects are dominated by the small particles. 
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terms have already been calculated (see, e.g., HH). The 
first steps are straight forward 



{iLr2Ti) = 




APPENDIX A: CALCULATION OF THE MIXED 
TERM 

We will only show the calculations for the mixed term 
{i^i2Ti)^ for D = 2. The calculations in D = 3 are very 
similar although slightly longer and the single species 



where we used the molecular chaos assumption to reduce 
the average over all possible pairs of colliding spheres 
to a sum of 2N1N2 times the average result of a single 
colliding pair. 

Now we introduce two partitions of unity 
(/ d^Rid^R2S{Ri - ri)(5(J?2 - ra)), i.e. 



UL12T1) = N2 



J dT J d^Rid^R25iR, - r,)S{R2 ~ r2)fNav^},t)i'J[f^'-^vl 
JdrM{v,},t) 



identifying the pair correlation function gi2{Ri2) /V^ = {6{Ri — ri)S{R2 — T2))t in this expression yields 

(i£;12Tl) 



N2 I d^Rid^R-i J R d\g,2{Ri2)fN{{v^}, t)i7[f^'-^vl 



y2 



/n,d2«,/Ar(K},t) 



Substituting R12 = Ri — R2 for Ri the other spatial integration is trivial as are all the velocity integrals in the 
denominator and those for j > 2 in the numerator: 



(«iLi2Ti) = —X2n 



nil 



2^Ti{t) 



1712 



2TlT2{t) 



d'^i?i2 / d^wid^U2exp 









2Ti 


exp 


2Ti _ 



Writing R12 in polar coordinates such that R12 ■ V12 = W12 cos (p, the radial integration is simply the application of 



the ^-function in 7[^^' and the step function constraints the angular integration. 



37r/2 
it/2 



d^tiid^W2'yi2 cos (f> exp 



rriivf 




m2vl ' 




exp 


[ 2T2it)\ 



where 



"N = X2n{ai + 02)% 



12 



mi 



2T:Ti{t) 



1712 



27rT2(t) 



According to the collision rules, the application of b'{2^ yields 



ib^i2^ - 1)^? = - — (1 + ei2)(n • ^i2)(n ■vi) + ^{l + eufin ■ Vi2f 
mi mf 



where fi = /ii2 is the reduced mass. Introducing the new average 

f3TT/2 



7r/2 



d^i;id^i;2Ui2 cos (j>A exp 







m2t;| " 


[ 2Ti(t)J 


exp 


[ 2T2{t)\ 
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what we have to calculate is 



(iiLiaTi) = + ei2) ((n • Vi2){fi ■ vi))^ + £^^(1 + £12)' ((n • vi2f)^ 



(Al) 



Let's consider the first term in eq. (|Aip . Substituting v = V12 for V2 and writing Vi in polar coordinates such that 
Vi ■ V ~ wit;cos7 one gets 

^37r/2 coo /.27r 



((1)12 ■ n)ivi ■ h))^ ^ d V 



d(j) 



7r/2 Jo JQ 



X exp 



1 m2Ti{t) + miT2{t) 



COS^ (, 


/)cos(7 — 








m2V^ 




7712^^1''^ 


exp 


[ 2T2(t)J 


exp 


7— COS 7 

[ 7^2 (t) 'J 



! Ti{t)T2{t) 

Invoking the addition theorem for cos (7 — (j)) the integration over (j) becomes trivial and the integration over 7 defines 
the associated Bcssel function Ii [x) . 



((■"12 ■n)ivi •n))2 = — — J d^v J dviv'^vlli{m2VVi/T2{t)) exp 



1 m2Ti{t) + miT2{t) 



2 Ti{t)T2{t) 

Integrals of the form j dxx^^^ In{oix) exp(— /3a::^) have closed solutions such that we get 





m2V^ 


exp 


[ 2T2{t)\ 



{{vri ■ n){vi ■ n))2 = - 



Stt 7712 

YT2{t) 



Ti{t)T2{t) 



vv exp 



mim2 



2 m2Ti{t) + miT2{t) 



_m2Ti{t) + miT2{t) _ 

We are left with a pair of Gaussian integrals. Calculating the second term in eq. (|A1[) involves essentially the same 
steps as shown above. 



APPENDIX B: SOLVING THE INTEGRAL 
EQUATIONS 



To be able to apply Banach's fix point iteration, we 
rearrange eqs. pO]) and (fT2|) and define operators 



-G / + rf V'^ + ^(1 + ef^T{r)dr H{a) 
MT]{a) := 



C/;^(a + r)y2M + IM (,2 _ 1) _ (1 + ,)2_^ 



dr 



and 



A2H '](«) := 



-C J^{a + r)V^ + ^(l + ef^u^-\r)dr 2uj^'{a) 
Ri 



Ri 



Cj^{a + r)^^J^ + ^ (,2_i)_(i + ,)2_^ 



dr 



with C = ?i\/6/(i?2 — Ri)- Now the solutions of the 
integral equations are the fix points of Ai and A2 , which 
we try to determine by iteration. This method worked 
well in the case of eq. ^TU\i . for ujo, however, convergence 



was not fully satisfactory. 

That is why we combined it with Newtons method. We 
define the function G whose root is to be determined by 



G[f]{a) = 2f{a) + C 



Aa + rf 



/(a) , fir) 



{e^^l)f{a) + il + ef 



dr 



Ri 
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and calculate its functional derivative. After discretiza- 
tion of the integrals we obtain a function G : R*'' M,^^ 
on which we can apply Newton's method. Newton's 



method requiring a sufficiently good starting approxima- 
tion, we chose as such the result of Banach's fixpoint 
iteration after about 300 iterations. 
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